simion.workbench_program()

adjustable simion_output_mode = 0 
adjustable pe_update_each_usec = 0.01 
adjustable _rfvolt = 2500 
adjustable _rfvolt2 = 1000 
adjustable _frequency = 10000
adjustable _ke = 90000

adjustable _length = 500

local last_pe_update = 0.0

local PLOT = simion.import(
  ((PLOT_LIBRARY or 'excel') == 'excel') and 'excellib.lua' or
  (PLOT_LIBRARY == 'gnuplot')            and '../gnuplot/gnuplotlib.lua' or
  error('invalid PLOT_LIBRARY: ' .. tostring(PLOT_LIBRARY))
)

local particle_count = 0
local dataset = {}

local function set_frequency(frequency)                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
	local omega = 2 * math.pi * frequency * 1E-6	
	return omega
end

function segment.fast_adjust()

	local omega = set_frequency(_frequency) 
	
	
	adj_elect01 = _rfvolt * sin(omega * ion_time_of_flight)
	adj_elect02 = _rfvolt2 * sin(omega * ion_time_of_flight + 1 * math.pi / 2)
	adj_elect03 = _rfvolt * sin(omega * ion_time_of_flight + 2 * math.pi / 2)
	adj_elect04 = _rfvolt2 * sin(omega * ion_time_of_flight + 3 * math.pi / 2)
	
	adj_elect05 = _rfvolt * sin(omega * ion_time_of_flight +  math.pi)
	adj_elect06 = _rfvolt2 * sin(omega * ion_time_of_flight + 1 * math.pi / 2 + math.pi)
	adj_elect07 = _rfvolt * sin(omega * ion_time_of_flight +  2 * math.pi / 2 + math.pi)
	adj_elect08 = _rfvolt2 * sin(omega * ion_time_of_flight + 3 * math.pi / 2 + math.pi)
end	

function segment.other_actions()
	if abs(ion_time_of_flight - last_pe_update) >= pe_update_each_usec then
	last_pe_update = ion_time_of_flight 
	sim_update_pe_surface = 1
	end

	if ion_splat == 0 then return end

	local x, y, z = ion_px_mm, ion_py_mm, ion_pz_mm
	local mass = ion_mass
	local vx, vy, vz = ion_vx_mm, ion_vy_mm, ion_vz_mm
	local t = ion_time_of_flight
	local charge = ion_charge
	local f = _frequency
	
	-- KE = 0.5 * mass * v^2 × 转换系数
	-- mass: amu, v: mm/μs, KE: eV
	-- 转换系数: (1.66054e-27 kg/amu) × (1e6 m²/s² per (mm/μs)²) / (1.60218e-19 J/eV) = 0.01036427
	local AMU_MMPS2_TO_EV = 0.01036427
	local v_squared = vx*vx + vy*vy + vz*vz
	local ke = 0.5 * mass * v_squared * AMU_MMPS2_TO_EV

	print(string.format("Splat at x=%.6f mm, y=%.6f mm, z=%.6f mm, mass=%.6f amu, vx=%.6f mm/us, vy=%.6f mm/us, vz=%.6f mm/us, KE=%.6f eV, t=%.6f us, charge=%.6f, f=%.6f", x, y, z, mass, vx, vy, vz, ke, t, charge, f))
	particle_count = particle_count + 1
    dataset[particle_count] = {x, y, z, mass, f, vx, vy, vz, ke, t, charge}
end

function segment.tstep_adjust()
	ion_time_step = min(ion_time_step, 0.001) 
end

function segment.terminate_run()
	-- local f = assert(io.open(_frequency .. ".csv", "w"))
	-- f:write("x,y,z,mass,f\n")
  	-- for _, row in ipairs(dataset) do
    	-- f:write(string.format("%.6f,%.6f,%.6f,%.6f,%.6f\n", row[1], row[2], row[3], row[4], row[5]))
  	-- end
  	-- f:close()	
    -- PLOT.plot {title='SIMION Phase Plot', xlabel='y', ylabel='yprime', dataset}
end